% Extended Merton model analysis of climate risk effects on municipal bond yields.
% Mike Schwert 10/6/21

clear all;

% set directory for figure output
% cd 'C:\Users\schwert\Dropbox (Personal)\GGLS_munis\Paper\Figures';

% colors for plots
navy = [0 0 (128/255)]; maroon = [(128/255) 0 0]; darkgreen = [0 (100/255) 0]; darkorange = [(255/256) (140/255) 0];

% set baseline parameters
global y t r a v k dy dv vol0;
options = optimset('TolFun', 1.0e-10, 'MaxIter', 100000, 'TolX',1.0e-10,'MaxFunEvals', 100000, 'Display', 'iter');

rf = 0.0268;        % risk-free rate
lmda = 0;           % jump intensity (Merton 1976)
r = rf + lmda;      % adjusted risk-free rate (Merton 1976)
y = rf + 0.0056;    % bond yield - 56 bps spread (robustness to tax-adjustment: 190 bps)
t = 7.5;            % typical duration (vs. 10-year maturity; robustness: 22.5 years)
v = 100;            % value of assets (unobservable for munis - set to arbitrary level)
k = 10;             % face value of debt (difficult to observe - vary with asset value)

% find the set of (volatility, leverage) combinations that match the yield
lev_vol = zeros(97,2); errors = zeros(97,1);
for i = 2:98
    k = i;
    lev_vol(i-1,1) = i;
    [vol, errors(i)] = fminsearch(@ggls_yield_opt, 0.20, options);
    lev_vol(i-1,2) = vol*100;
end

figure1 = figure('WindowState','maximized','OuterPosition',[-1.003515625 0.0284722222222222 1.00703125 0.977777777777778]);
axes1 = axes('Parent',figure1); hold(axes1,'on');
plot1 = plot(lev_vol(:,1), lev_vol(:,2), 'LineStyle', '-', 'Color', navy, 'LineWidth', 2);
ylabel('Asset Volatility (%)','FontName','Times New Roman','FontSize',32);
xlabel('Leverage Ratio (%)','FontName','Times New Roman','FontSize',32);
box(axes1,'on'); set(axes1,'FontName','Times New Roman','FontSize',32);
saveas(gcf, 'vol_lev_calc_wide.png');

% select a few levels of leverage and vary the shock to V
lev_yields = zeros(260,5);
for i = 0:259
    lev_yields(i+1,1) = i/10;
    lev_yields(i+1,2) = 100*(ggls_yield_calc(t, r, (1 - i/1000)*v, lev_vol(10,1), lev_vol(10,2)/100) - ggls_yield_calc(t, r, v, lev_vol(10,1), lev_vol(10,2)/100));
    lev_yields(i+1,3) = 100*(ggls_yield_calc(t, r, (1 - i/1000)*v, lev_vol(30,1), lev_vol(30,2)/100) - ggls_yield_calc(t, r, v, lev_vol(30,1), lev_vol(30,2)/100));
    lev_yields(i+1,4) = 100*(ggls_yield_calc(t, r, (1 - i/1000)*v, lev_vol(50,1), lev_vol(50,2)/100) - ggls_yield_calc(t, r, v, lev_vol(50,1), lev_vol(50,2)/100));
    lev_yields(i+1,5) = 100*(ggls_yield_calc(t, r, (1 - i/1000)*v, lev_vol(70,1), lev_vol(70,2)/100) - ggls_yield_calc(t, r, v, lev_vol(70,1), lev_vol(70,2)/100));
end

figure1 = figure('WindowState','maximized','OuterPosition',[-1.003515625 0.0284722222222222 1.00703125 0.977777777777778]);
axes1 = axes('Parent',figure1); hold(axes1,'on');
box(axes1,'on'); xlim(axes1,[0 25]); set(axes1, 'FontName', 'Times New Roman', 'FontSize', 32, 'XTick', [0 5 10 15 20 25]);
plot1 = plot(lev_yields(:,1), lev_yields(:,2), 'LineStyle', '-', 'Color', navy, 'LineWidth', 2);
plot2 = plot(lev_yields(:,1), lev_yields(:,3), 'LineStyle', '--', 'Color', maroon, 'LineWidth', 2);
plot3 = plot(lev_yields(:,1), lev_yields(:,4), 'LineStyle', '-.', 'Color', darkgreen, 'LineWidth', 2);
plot4 = plot(lev_yields(:,1), lev_yields(:,5), 'LineStyle', ':', 'Color', darkorange, 'LineWidth', 2);
ylabel('Change in Yield (%)', 'FontName', 'Times New Roman', 'FontSize', 32);
xlabel('Decrease in PV(Cash Flows) (%)', 'FontName', 'Times New Roman', 'FontSize', 32);
legend1 = legend('K/V = 0.10', 'K/V = 0.30', 'K/V = 0.50', 'K/V = 0.70');
set(legend1, 'Position', [0.157779945886384 0.713573812898783 0.110058595780283 0.171971454590927], 'EdgeColor', [1 1 1]);
saveas(gcf, 'dy_dlev_baseline_wide.png');

% check on the effect of an increase in volatility
vol_yields = zeros(260,5);
for i = 0:259
    vol_yields(i+1,1) = i/10;
    vol_yields(i+1,2) = 100*(ggls_yield_calc(t, r, v, lev_vol(10,1), (1 + i/1000)*lev_vol(10,2)/100) - ggls_yield_calc(t, r, v, lev_vol(10,1), lev_vol(10,2)/100));
    vol_yields(i+1,3) = 100*(ggls_yield_calc(t, r, v, lev_vol(30,1), (1 + i/1000)*lev_vol(30,2)/100) - ggls_yield_calc(t, r, v, lev_vol(30,1), lev_vol(30,2)/100));
    vol_yields(i+1,4) = 100*(ggls_yield_calc(t, r, v, lev_vol(50,1), (1 + i/1000)*lev_vol(50,2)/100) - ggls_yield_calc(t, r, v, lev_vol(50,1), lev_vol(50,2)/100));
    vol_yields(i+1,5) = 100*(ggls_yield_calc(t, r, v, lev_vol(70,1), (1 + i/1000)*lev_vol(70,2)/100) - ggls_yield_calc(t, r, v, lev_vol(70,1), lev_vol(70,2)/100));
end

figure1 = figure('WindowState','maximized','OuterPosition',[-1.003515625 0.0284722222222222 1.00703125 0.977777777777778]);
axes1 = axes('Parent',figure1); hold(axes1,'on');
box(axes1,'on'); xlim(axes1,[0 25]); set(axes1, 'FontName', 'Times New Roman', 'FontSize', 32, 'XTick', [0 5 10 15 20 25]);
plot1 = plot(vol_yields(:,1), vol_yields(:,2), 'LineStyle', '-', 'Color', navy, 'LineWidth', 2);
plot2 = plot(vol_yields(:,1), vol_yields(:,3), 'LineStyle', '--', 'Color', maroon, 'LineWidth', 2);
plot3 = plot(vol_yields(:,1), vol_yields(:,4), 'LineStyle', '-.', 'Color', darkgreen, 'LineWidth', 2);
plot4 = plot(vol_yields(:,1), vol_yields(:,5), 'LineStyle', ':', 'Color', darkorange, 'LineWidth', 2);
ylabel('Change in Yield (%)','FontName','Times New Roman','FontSize',32);
xlabel('Increase in Volatility of PV(Cash Flows) (%)','FontName','Times New Roman','FontSize',32);
legend1 = legend('K/V = 0.10', 'K/V = 0.30', 'K/V = 0.50', 'K/V = 0.70');
set(legend1, 'Position', [0.157779945886384 0.713573812898783 0.110058595780283 0.171971454590927], 'EdgeColor', [1 1 1]);
saveas(gcf, 'dy_dvol_baseline_wide.png');

% find the combination of V and and volatility that matches the change in yield implied by our regressions
dy = 0.00053;       % estimated effect of SLR exposure - 5.3 bps
y = rf + 0.0063;    % bond yield - 63 bps spread
t = 7.5;            % typical duration (vs. 10-year maturity; robustness: 22.5 years)

dv_dvol_ggls = zeros(8,5); errors_ggls = zeros(8,4);

for i = -5:2
    dv = i/100; dv_dvol_ggls(i+6, 1) = dv*100;
    
    k = 10; vol0 = lev_vol(10,2)/100;
    [dvol, errors_ggls(i+6,1)] = fminsearch(@ggls_dy_dvol_opt, 0.01, options);
    dv_dvol_ggls(i+6,2) = dvol*100;
    
    k = 30; vol0 = lev_vol(30,2)/100;
    [dvol, errors_ggls(i+6,2)] = fminsearch(@ggls_dy_dvol_opt, 0.01, options);
    dv_dvol_ggls(i+6,3) = dvol*100;
    
    k = 50; vol0 = lev_vol(50,2)/100;
    [dvol, errors_ggls(i+6,3)] = fminsearch(@ggls_dy_dvol_opt, 0.01, options);
    dv_dvol_ggls(i+6,4) = dvol*100;
    
    k = 70; vol0 = lev_vol(70,2)/100;
    [dvol, errors_ggls(i+6,4)] = fminsearch(@ggls_dy_dvol_opt, 0.01, options);
    dv_dvol_ggls(i+6,5) = dvol*100;
end

figure1 = figure('WindowState','maximized','OuterPosition',[-1.003515625 0.0284722222222222 1.00703125 0.977777777777778]);
axes1 = axes('Parent',figure1); hold(axes1,'on');
plot1 = plot(dv_dvol_ggls(:,1), dv_dvol_ggls(:,2), 'LineStyle', '-', 'Color', navy, 'LineWidth', 2);
plot2 = plot(dv_dvol_ggls(:,1), dv_dvol_ggls(:,3), 'LineStyle', '--', 'Color', maroon, 'LineWidth', 2);
plot3 = plot(dv_dvol_ggls(:,1), dv_dvol_ggls(:,4), 'LineStyle', '-.', 'Color', darkgreen, 'LineWidth', 2);
plot4 = plot(dv_dvol_ggls(:,1), dv_dvol_ggls(:,5), 'LineStyle', ':', 'Color', darkorange, 'LineWidth', 2);
ylabel('Change in Volatility of PV(CF) (%)','FontName','Times New Roman','FontSize',32);
xlabel('Change in PV(Cash Flows) (%)','FontName','Times New Roman','FontSize',32);
box(axes1,'on'); set(axes1,'FontName','Times New Roman','FontSize',32);
legend1 = legend('K/V = 0.10', 'K/V = 0.30', 'K/V = 0.50', 'K/V = 0.70');
set(legend1, 'Position', [0.157779945886384 0.713573812898783 0.110058595780283 0.171971454590927], 'EdgeColor', [1 1 1]);
saveas(gcf, 'dv_dvol_ggls_wide.png');

% find the change in volatility implied by a 1.4% reduction in V
dv = -0.014; dvol_v14 = zeros(4,1);
k = 10; vol0 = lev_vol(10,2)/100;
[dvol_v14(1), e] = fminsearch(@ggls_dy_dvol_opt, 0.01, options);
k = 30; vol0 = lev_vol(30,2)/100;
[dvol_v14(2), e] = fminsearch(@ggls_dy_dvol_opt, 0.01, options);
k = 50; vol0 = lev_vol(50,2)/100;
[dvol_v14(3), e] = fminsearch(@ggls_dy_dvol_opt, 0.01, options);
k = 70; vol0 = lev_vol(70,2)/100;
[dvol_v14(4), e] = fminsearch(@ggls_dy_dvol_opt, 0.01, options);

% do the same thing for the main result in Painter (2020)
dy = 0.00234;           % estimated effect of climate exposure - 23.4 bps
rf = 0.040;             % higher risk-free rate at longer maturity
y = rf + 0.0070;        % bond yield - 63 bps spread
r = rf + lmda;          % adjusted risk-free rate (Merton 1976)
t = 22.5;               % typical duration of 30-year bond

dv_dvol_paint = zeros(36,5); errors_paint = zeros(36,4);

for i = -35:0
    dv = i/100; dv_dvol_paint(i+36, 1) = dv*100;
    
    k = 10; vol0 = lev_vol(10,2)/100;
    [dvol, errors_paint(i+36,1)] = fminsearch(@ggls_dy_dvol_opt, 0.01, options);
    dv_dvol_paint(i+36,2) = dvol*100;
    
    k = 30; vol0 = lev_vol(30,2)/100;
    [dvol, errors_paint(i+36,2)] = fminsearch(@ggls_dy_dvol_opt, 0.01, options);
    dv_dvol_paint(i+36,3) = dvol*100;
    
    k = 50; vol0 = lev_vol(50,2)/100;
    [dvol, errors_paint(i+36,3)] = fminsearch(@ggls_dy_dvol_opt, 0.01, options);
    dv_dvol_paint(i+36,4) = dvol*100;
    
    k = 70; vol0 = lev_vol(70,2)/100;
    [dvol, errors_paint(i+36,4)] = fminsearch(@ggls_dy_dvol_opt, 0.01, options);
    dv_dvol_paint(i+36,5) = dvol*100;
end

figure1 = figure('WindowState','maximized','OuterPosition',[-1.003515625 0.0284722222222222 1.00703125 0.977777777777778]);
axes1 = axes('Parent',figure1); hold(axes1,'on');
plot1 = plot(dv_dvol_paint(:,1), dv_dvol_paint(:,2), 'LineStyle', '-', 'Color', navy, 'LineWidth', 2);
plot2 = plot(dv_dvol_paint(:,1), dv_dvol_paint(:,3), 'LineStyle', '--', 'Color', maroon, 'LineWidth', 2);
plot3 = plot(dv_dvol_paint(:,1), dv_dvol_paint(:,4), 'LineStyle', '-.', 'Color', darkgreen, 'LineWidth', 2);
plot4 = plot(dv_dvol_paint(:,1), dv_dvol_paint(:,5), 'LineStyle', ':', 'Color', darkorange, 'LineWidth', 2);
ylabel('Change in Volatility of PV(CF) (%)','FontName','Times New Roman','FontSize',32);
xlabel('Change in PV(Cash Flows) (%)','FontName','Times New Roman','FontSize',32);
box(axes1,'on'); set(axes1,'FontName','Times New Roman','FontSize',32);
legend1 = legend('K/V = 0.10', 'K/V = 0.30', 'K/V = 0.50', 'K/V = 0.70');
set(legend1, 'Position', [0.157779945886384 0.713573812898783 0.110058595780283 0.171971454590927], 'EdgeColor', [1 1 1]);
saveas(gcf, 'dv_dvol_painter_wide.png');

% yield sensitivity to shocks when bankruptcy costs differ across states
t = 7.5; rf = 0.0268; lmda = 0; r = rf + lmda; y = rf + 0.0056; v = 100; k = 40;

bc_vol = zeros(51,2); errors = zeros(51,1);
for i = 1:51
    a = (i-1)/100;
    bc_vol(i,1) = a;
    [vol, errors(i)] = fminsearch(@ggls_yield_opt_bc, 0.20, options);
    bc_vol(i,2) = vol*100;
end

lev_yields_bc = zeros(260,7);
for i = 0:259
    lev_yields_bc(i+1,1) = i/10;
    lev_yields_bc(i+1,2) = 100*(ggls_yield_calc_bc(t, r, bc_vol(1,1), (1 - i/1000)*v, k, bc_vol(1,2)/100) - ggls_yield_calc_bc(t, r, bc_vol(1,1), v, k, bc_vol(1,2)/100));
    lev_yields_bc(i+1,3) = 100*(ggls_yield_calc_bc(t, r, bc_vol(11,1), (1 - i/1000)*v, k, bc_vol(11,2)/100) - ggls_yield_calc_bc(t, r, bc_vol(11,1), v, k, bc_vol(11,2)/100));
    lev_yields_bc(i+1,4) = 100*(ggls_yield_calc_bc(t, r, bc_vol(21,1), (1 - i/1000)*v, k, bc_vol(21,2)/100) - ggls_yield_calc_bc(t, r, bc_vol(21,1), v, k, bc_vol(21,2)/100));
    lev_yields_bc(i+1,5) = 100*(ggls_yield_calc_bc(t, r, bc_vol(31,1), (1 - i/1000)*v, k, bc_vol(31,2)/100) - ggls_yield_calc_bc(t, r, bc_vol(31,1), v, k, bc_vol(31,2)/100));
    lev_yields_bc(i+1,6) = 100*(ggls_yield_calc_bc(t, r, bc_vol(41,1), (1 - i/1000)*v, k, bc_vol(41,2)/100) - ggls_yield_calc_bc(t, r, bc_vol(41,1), v, k, bc_vol(41,2)/100));
    lev_yields_bc(i+1,7) = 100*(ggls_yield_calc_bc(t, r, bc_vol(51,1), (1 - i/1000)*v, k, bc_vol(51,2)/100) - ggls_yield_calc_bc(t, r, bc_vol(51,1), v, k, bc_vol(51,2)/100));
end

figure1 = figure('WindowState','maximized','OuterPosition',[-1.003515625 0.0284722222222222 1.00703125 0.977777777777778]);
axes1 = axes('Parent',figure1); hold(axes1,'on');
box(axes1,'on'); xlim(axes1,[0 25]); set(axes1, 'FontName', 'Times New Roman', 'FontSize', 32, 'XTick', [0 5 10 15 20 25]);
plot1 = plot(lev_yields_bc(:,1), lev_yields_bc(:,2), 'LineStyle', '-', 'Color', navy, 'LineWidth', 2);
plot2 = plot(lev_yields_bc(:,1), lev_yields_bc(:,3), 'LineStyle', '--', 'Color', maroon, 'LineWidth', 2);
plot3 = plot(lev_yields_bc(:,1), lev_yields_bc(:,4), 'LineStyle', '-.', 'Color', darkgreen, 'LineWidth', 2);
plot4 = plot(lev_yields_bc(:,1), lev_yields_bc(:,5), 'LineStyle', ':', 'Color', darkorange, 'LineWidth', 2);
ylabel('Change in Yield (%)','FontName','Times New Roman','FontSize',32);
xlabel('Decrease in PV(Cash Flows) (%)','FontName','Times New Roman','FontSize',32);
box(axes1,'on'); set(axes1,'FontName','Times New Roman','FontSize',32);
legend1 = legend('a = 0', 'a = 0.10', 'a = 0.20', 'a = 0.30');
set(legend1, 'Position', [0.157779945886384 0.713573812898783 0.110058595780283 0.171971454590927], 'EdgeColor', [1 1 1]);
saveas(gcf, 'dy_dlev_distress_wide.png');

vol_yields_bc = zeros(260,7);
for i = 0:259
    vol_yields_bc(i+1,1) = i/10;
    vol_yields_bc(i+1,2) = 100*(ggls_yield_calc_bc(t, r, bc_vol(1,1), v, k, (1 + i/1000)*bc_vol(1,2)/100) - ggls_yield_calc_bc(t, r, bc_vol(1,1), v, k, bc_vol(1,2)/100));
    vol_yields_bc(i+1,3) = 100*(ggls_yield_calc_bc(t, r, bc_vol(11,1), v, k, (1 + i/1000)*bc_vol(11,2)/100) - ggls_yield_calc_bc(t, r, bc_vol(11,1), v, k, bc_vol(11,2)/100));
    vol_yields_bc(i+1,4) = 100*(ggls_yield_calc_bc(t, r, bc_vol(21,1), v, k, (1 + i/1000)*bc_vol(21,2)/100) - ggls_yield_calc_bc(t, r, bc_vol(21,1), v, k, bc_vol(21,2)/100));
    vol_yields_bc(i+1,5) = 100*(ggls_yield_calc_bc(t, r, bc_vol(31,1), v, k, (1 + i/1000)*bc_vol(31,2)/100) - ggls_yield_calc_bc(t, r, bc_vol(31,1), v, k, bc_vol(31,2)/100));
    vol_yields_bc(i+1,6) = 100*(ggls_yield_calc_bc(t, r, bc_vol(41,1), v, k, (1 + i/1000)*bc_vol(41,2)/100) - ggls_yield_calc_bc(t, r, bc_vol(41,1), v, k, bc_vol(41,2)/100));
    vol_yields_bc(i+1,7) = 100*(ggls_yield_calc_bc(t, r, bc_vol(51,1), v, k, (1 + i/1000)*bc_vol(51,2)/100) - ggls_yield_calc_bc(t, r, bc_vol(51,1), v, k, bc_vol(51,2)/100));
end

figure1 = figure('WindowState','maximized','OuterPosition',[-1.003515625 0.0284722222222222 1.00703125 0.977777777777778]);
axes1 = axes('Parent',figure1); hold(axes1,'on');
box(axes1,'on'); xlim(axes1,[0 25]); set(axes1, 'FontName', 'Times New Roman', 'FontSize', 32, 'XTick', [0 5 10 15 20 25]);
plot1 = plot(vol_yields_bc(:,1), vol_yields_bc(:,2), 'LineStyle', '-', 'Color', navy, 'LineWidth', 2);
plot2 = plot(vol_yields_bc(:,1), vol_yields_bc(:,3), 'LineStyle', '--', 'Color', maroon, 'LineWidth', 2);
plot3 = plot(vol_yields_bc(:,1), vol_yields_bc(:,4), 'LineStyle', '-.', 'Color', darkgreen, 'LineWidth', 2);
plot4 = plot(vol_yields_bc(:,1), vol_yields_bc(:,5), 'LineStyle', ':', 'Color', darkorange, 'LineWidth', 2);
ylabel('Change in Yield (%)','FontName','Times New Roman','FontSize',32);
xlabel('Increase in Volatility of PV(Cash Flows) (%)','FontName','Times New Roman','FontSize',32);
box(axes1,'on'); set(axes1,'FontName','Times New Roman','FontSize',32);
legend1 = legend('a = 0', 'a = 0.10', 'a = 0.20', 'a = 0.30');
set(legend1, 'Position', [0.157779945886384 0.713573812898783 0.110058595780283 0.171971454590927], 'EdgeColor', [1 1 1]);
saveas(gcf, 'dy_dvol_distress_wide.png');

% yield sensitivity to shocks when the issuer can have senior bank debt
t = 7.5; rf = 0.0268; lmda = 0; r = rf + lmda; y = rf + 0.0056; v = 100; k = 40;

jr_vol = zeros(76,2); errors = zeros(76,1);
for i = 1:76
    s_k = (i-1)/100;
    jr_vol(i,1) = s_k;
    [vol, errors(i)] = fminsearch(@ggls_yield_opt_jr, 0.20, options);
    jr_vol(i,2) = vol*100;
end

lev_yields_jr = zeros(260,5);
for i = 0:259
    lev_yields_jr(i+1,1) = i/10;
    lev_yields_jr(i+1,2) = 100*(ggls_yield_calc_jr(t, r, (1 - i/1000)*v, k, jr_vol(1,1), jr_vol(1,2)/100, 0) - ggls_yield_calc_jr(t, r, v, k, jr_vol(1,1), jr_vol(1,2)/100, 0));
    lev_yields_jr(i+1,3) = 100*(ggls_yield_calc_jr(t, r, (1 - i/1000)*v, k, jr_vol(26,1), jr_vol(26,2)/100, 0) - ggls_yield_calc_jr(t, r, v, k, jr_vol(26,1), jr_vol(26,2)/100, 0));
    lev_yields_jr(i+1,4) = 100*(ggls_yield_calc_jr(t, r, (1 - i/1000)*v, k, jr_vol(51,1), jr_vol(51,2)/100, 0) - ggls_yield_calc_jr(t, r, v, k, jr_vol(51,1), jr_vol(51,2)/100, 0));
    lev_yields_jr(i+1,5) = 100*(ggls_yield_calc_jr(t, r, (1 - i/1000)*v, k, jr_vol(76,1), jr_vol(76,2)/100, 0) - ggls_yield_calc_jr(t, r, v, k, jr_vol(76,1), jr_vol(76,2)/100, 0));
end

figure1 = figure('WindowState','maximized','OuterPosition',[-1.003515625 0.0284722222222222 1.00703125 0.977777777777778]);
axes1 = axes('Parent',figure1); hold(axes1,'on');
box(axes1,'on'); xlim(axes1,[0 25]); set(axes1, 'FontName', 'Times New Roman', 'FontSize', 32, 'XTick', [0 5 10 15 20 25]);
plot1 = plot(lev_yields_jr(:,1), lev_yields_jr(:,2), 'LineStyle', '-', 'Color', navy, 'LineWidth', 2);
plot2 = plot(lev_yields_jr(:,1), lev_yields_jr(:,3), 'LineStyle', '--', 'Color', maroon, 'LineWidth', 2);
plot3 = plot(lev_yields_jr(:,1), lev_yields_jr(:,4), 'LineStyle', '-.', 'Color', darkgreen, 'LineWidth', 2);
ylabel('Change in Yield (%)','FontName','Times New Roman','FontSize',32);
xlabel('Decrease in PV(Cash Flows) (%)','FontName','Times New Roman','FontSize',32);
box(axes1,'on'); set(axes1,'FontName','Times New Roman','FontSize',32);
legend1 = legend('K_S/K = 0', 'K_S/K = 0.25', 'K_S/K = 0.50');
set(legend1, 'Position', [0.157779945886384 0.713573812898783 0.110058595780283 0.171971454590927], 'EdgeColor', [1 1 1]);
saveas(gcf, 'dy_dlev_tiered_wide.png');

vol_yields_jr = zeros(260,5);
for i = 0:259
    vol_yields_jr(i+1,1) = i/10;
    vol_yields_jr(i+1,2) = 100*(ggls_yield_calc_jr(t, r, v, k, jr_vol(1,1), (1 + i/1000)*jr_vol(1,2)/100, 0) - ggls_yield_calc_jr(t, r, v, k, jr_vol(1,1), jr_vol(1,2)/100, 0));
    vol_yields_jr(i+1,3) = 100*(ggls_yield_calc_jr(t, r, v, k, jr_vol(11,1), (1 + i/1000)*jr_vol(11,2)/100, 0) - ggls_yield_calc_jr(t, r, v, k, jr_vol(11,1), jr_vol(11,2)/100, 0));
    vol_yields_jr(i+1,4) = 100*(ggls_yield_calc_jr(t, r, v, k, jr_vol(21,1), (1 + i/1000)*jr_vol(21,2)/100, 0) - ggls_yield_calc_jr(t, r, v, k, jr_vol(21,1), jr_vol(21,2)/100, 0));
    vol_yields_jr(i+1,5) = 100*(ggls_yield_calc_jr(t, r, v, k, jr_vol(31,1), (1 + i/1000)*jr_vol(31,2)/100, 0) - ggls_yield_calc_jr(t, r, v, k, jr_vol(31,1), jr_vol(31,2)/100, 0));
end

figure1 = figure('WindowState','maximized','OuterPosition',[-1.003515625 0.0284722222222222 1.00703125 0.977777777777778]);
axes1 = axes('Parent',figure1); hold(axes1,'on');
box(axes1,'on'); xlim(axes1,[0 25]); set(axes1, 'FontName', 'Times New Roman', 'FontSize', 32, 'XTick', [0 5 10 15 20 25]);
plot1 = plot(vol_yields_jr(:,1), vol_yields_jr(:,2), 'LineStyle', '-', 'Color', navy, 'LineWidth', 2);
plot2 = plot(vol_yields_jr(:,1), vol_yields_jr(:,3), 'LineStyle', '--', 'Color', maroon, 'LineWidth', 2);
plot3 = plot(vol_yields_jr(:,1), vol_yields_jr(:,4), 'LineStyle', '-.', 'Color', darkgreen, 'LineWidth', 2);
ylabel('Change in Yield (%)','FontName','Times New Roman','FontSize',32);
xlabel('Increase in Volatility of PV(Cash Flows) (%)','FontName','Times New Roman','FontSize',32);
box(axes1,'on'); set(axes1,'FontName','Times New Roman','FontSize',32);
legend1 = legend('K_S/K = 0', 'K_S/K = 0.25', 'K_S/K = 0.50');
set(legend1, 'Position', [0.157779945886384 0.713573812898783 0.110058595780283 0.171971454590927], 'EdgeColor', [1 1 1]);
saveas(gcf, 'dy_dvol_tiered_wide.png');


%%% tax-adjusted analysis %%%

% tax-adjusted version of the main analysis, following Schwert (2017)
    % adjust yield upwards and estimate parameters for taxable yields
    % analyze effects of shocks and adjust yields back to tax-exempt

tax_factor = 1/((1-0.35)*(1-0.05)); % average state tax rate is 5%

rf = 0.0310;                % risk-free rate from swap curve
lmda = 0;                   % jump intensity (Merton 1976)
r = rf + lmda;              % adjusted risk-free rate (Merton 1976)
y_e = 0.0324;               % tax-exempt bond yield
y = y_e*tax_factor;         % tax-adjusted yield
t = 7.5;                    % typical duration (vs. 10-year maturity; robustness: 22.5 years)
v = 100;                    % value of assets (unobservable for munis - set to arbitrary level)
k = 10;                     % face value of debt (difficult to observe - vary with asset value)

% find the set of (volatility, leverage) combinations that match the yield
lev_vol = zeros(97,2); errors = zeros(97,1);
for i = 2:98
    k = i;
    lev_vol(i-1,1) = i;
    [vol, errors(i)] = fminsearch(@ggls_yield_opt, 0.20, options);
    lev_vol(i-1,2) = vol*100;
end

figure1 = figure('WindowState','maximized','OuterPosition',[-1.003515625 0.0284722222222222 1.00703125 0.977777777777778]);
axes1 = axes('Parent',figure1); hold(axes1,'on');
plot1 = plot(lev_vol(:,1), lev_vol(:,2), 'LineStyle', '-', 'Color', navy, 'LineWidth', 2);
ylabel('Asset Volatility (%)','FontName','Times New Roman','FontSize',32);
xlabel('Leverage Ratio (%)','FontName','Times New Roman','FontSize',32);
box(axes1,'on'); set(axes1,'FontName','Times New Roman','FontSize',32);
saveas(gcf, 'vol_lev_calc_taxadj_wide.png');

% select a few levels of leverage and vary the shock to V
lev_yields = zeros(26,5);
for i = 0:25
    lev_yields(i+1,1) = i;
    lev_yields(i+1,2) = 100*(ggls_yield_calc(t, r, (1 - i/100)*v, lev_vol(10,1), lev_vol(10,2)/100)/tax_factor - y_e);
    lev_yields(i+1,3) = 100*(ggls_yield_calc(t, r, (1 - i/100)*v, lev_vol(30,1), lev_vol(30,2)/100)/tax_factor - y_e);
    lev_yields(i+1,4) = 100*(ggls_yield_calc(t, r, (1 - i/100)*v, lev_vol(50,1), lev_vol(50,2)/100)/tax_factor - y_e);
    lev_yields(i+1,5) = 100*(ggls_yield_calc(t, r, (1 - i/100)*v, lev_vol(70,1), lev_vol(70,2)/100)/tax_factor - y_e);
end

figure1 = figure('WindowState','maximized','OuterPosition',[-1.003515625 0.0284722222222222 1.00703125 0.977777777777778]);
axes1 = axes('Parent',figure1); hold(axes1,'on');
plot1 = plot(lev_yields(:,1), lev_yields(:,2), 'LineStyle', '-', 'Color', navy, 'LineWidth', 2);
plot2 = plot(lev_yields(:,1), lev_yields(:,3), 'LineStyle', '--', 'Color', maroon, 'LineWidth', 2);
plot3 = plot(lev_yields(:,1), lev_yields(:,4), 'LineStyle', '-.', 'Color', darkgreen, 'LineWidth', 2);
plot4 = plot(lev_yields(:,1), lev_yields(:,5), 'LineStyle', ':', 'Color', darkorange, 'LineWidth', 2);
ylabel('Change in Yield (%)','FontName','Times New Roman','FontSize',32);
xlabel('Decrease in PV(Cash Flows) (%)','FontName','Times New Roman','FontSize',32);
box(axes1,'on'); set(axes1,'FontName','Times New Roman','FontSize',32);
legend1 = legend('K/V = 0.10', 'K/V = 0.30', 'K/V = 0.50', 'K/V = 0.70');
set(legend1, 'Position', [0.157779945886384 0.713573812898783 0.110058595780283 0.171971454590927], 'EdgeColor', [1 1 1]);
saveas(gcf, 'dy_dlev_taxadj_wide.png');

% check on the effect of an increase in volatility
vol_yields = zeros(26,5);
for i = 0:25
    vol_yields(i+1,1) = i;
    vol_yields(i+1,2) = 100*(ggls_yield_calc(t, r, v, lev_vol(10,1), (1 + i/100)*lev_vol(10,2)/100)/tax_factor - y_e);
    vol_yields(i+1,3) = 100*(ggls_yield_calc(t, r, v, lev_vol(30,1), (1 + i/100)*lev_vol(30,2)/100)/tax_factor - y_e);
    vol_yields(i+1,4) = 100*(ggls_yield_calc(t, r, v, lev_vol(50,1), (1 + i/100)*lev_vol(50,2)/100)/tax_factor - y_e);
    vol_yields(i+1,5) = 100*(ggls_yield_calc(t, r, v, lev_vol(70,1), (1 + i/100)*lev_vol(70,2)/100)/tax_factor - y_e);
end

figure1 = figure('WindowState','maximized','OuterPosition',[-1.003515625 0.0284722222222222 1.00703125 0.977777777777778]);
axes1 = axes('Parent',figure1); hold(axes1,'on');
plot1 = plot(vol_yields(:,1), vol_yields(:,2), 'LineStyle', '-', 'Color', navy, 'LineWidth', 2);
plot2 = plot(vol_yields(:,1), vol_yields(:,3), 'LineStyle', '--', 'Color', maroon, 'LineWidth', 2);
plot3 = plot(vol_yields(:,1), vol_yields(:,4), 'LineStyle', '-.', 'Color', darkgreen, 'LineWidth', 2);
plot4 = plot(vol_yields(:,1), vol_yields(:,5), 'LineStyle', ':', 'Color', darkorange, 'LineWidth', 2);
ylabel('Change in Yield (%)','FontName','Times New Roman','FontSize',32);
xlabel('Increase in Volatility of PV(Cash Flows) (%)','FontName','Times New Roman','FontSize',32);
box(axes1,'on'); set(axes1,'FontName','Times New Roman','FontSize',32);
legend1 = legend('K/V = 0.10', 'K/V = 0.30', 'K/V = 0.50', 'K/V = 0.70');
set(legend1, 'Position', [0.157779945886384 0.713573812898783 0.110058595780283 0.171971454590927], 'EdgeColor', [1 1 1]);
saveas(gcf, 'dy_dvol_taxadj_wide.png');




